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Estimating the Probability of a Diffusing Target Encountering a 

Stationary Sensor 

by 

James N. Eagle 

Department of Operations Research 
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Monterey, CR 93943 

In this report two expressions are given for the probability of a 
diffusing target avoiding detection to time t by a sensor fixed at the 
center of a square region fl. The target is constrained to always remain 
within R. The first expression results from an approximation to the 
exact solution of the diffusion equation, and the second from 
experimentation with a Monte Carlo simulation of the diffusion process. 

The sensor considered is a definite range law dr "cookie cutter" 
detector. For such a sensor, there is a specified range, R, beyond which 
detection is impossible and inside which detection is certain. 

1. Background 

This work was begun with the intention of developing for the Naval 
UUar College, Newport, RI a simple expression for the probability of a 
diffusing target avoiding detection by a sensor conducting a moving, 
systematic search of the area R. It was soon realized that the special 
case of a stationary sensor had to be first addressed, and that this 
simpler problem was indeed nontrivial. 

The importance of the stationary sensor problem goes beyond being 
Limiting case of the moving sensor problem. It may be used to provide 
estimates of the rate at which randomly moving targets will encounter 



stationary objects with extended fields of influence, such as fixed 
acoustic sensors, sonobuoys, or possibly mines. 

. This stationary sensor problem might appear deceptively simple at 
first glance. Rfter all, it could be argued, this problem is equivalent to 
that of a searcher with a cookie cutter sensor of range R conducting a 
random search fora stationary target. Rnd Kooprnanl 1 9^6] and [1980] 
argued that the probability of nondetection to time t is exp(-2Rvt/R), 
where v is the speed of the random search. The problem is, of course, 
what to use for v when the searcher's path is a diffusion. One of the 
results of this work is an expression for the equivalent speed of such a 
"diffusion search" of a square area. 

The initial, experimental results for the problem here addressed were 
obtained by Sislioglu[1981]. Sislioglu conducted Monte Carlo simulations 
with different target diffusion constants D, area sizes R, and sensor 
detection ranges R. He observed that when the initial distribution of the 
target was uniform over R, the probability of nondetection to time t, 
PND(t), was given approximately by 

(1 - 7TR 2 /R) exp(-247 RDt/R 15 ). (1) 

In this report, some analytical support for Sislioglu's results is offered. 
Rlso a slight modification of (1) is suggested which agrees more closely 
with theory and experimental results when the area of the detection 
disk approaches the area of the search region R. , 

2. The Diffusion Equation 

The probability density p(t) of a particle undergoing diffusion in any 
coordinate system must satisfy the diffusion equation 

(D/2] V 2 p = dp/at. 



( 2 ) 



where D is the diffusion constant, and V 2 is the Laplacian operator. In 
Cartesian and polar coordinates, respectively, (2) becomes 

(D/2) (3 2 p/3x 2 + 3 2 p/3y 2 ) = 3p/3t, arid (3) 

(D/2) (3 2 p/3r 2 + (l/r 2 )0 2 p/36 2 ) + (1/r)(3p/3r)) = 3p/3t. (4) 

To find a unique function p(x,y,t) or p(r,8,t) satisfying these partial 
differential equations, spatial and temporal boundary conditions must be 
specified. Defining R to be a square with sides of length L in the first 
quadrant, the boundary conditions in Cartesian coordinates are 



3p/3x|x=0 = 3p/3x| x= L = 0, (5) 

3p/9y |y=0 = 3p/9y|y=L = 0 , (6) 

p(x,y,t) = 0, ((x - L/2) 2 +(y - L/2) 2 ) 1/2 « R, and (7) 

p(x,y,0) = 1/R, ((x - L/2) 2 +(y - L/2) 2 ) 1/2 > R. (8) 



Equations (5) and (6) ensure that none of the target's probability mass 
"escapes" from R. That is, the boundaries of R are reflecting. Equation 
(7) requires that p(x,y,t) be 0 on the detection disk. Rrid (8) ensures that 
the initial distribution of the target over the search region is uniform and 
integrates to (R-rrR 2 )/R (the probability that the target is not detected 
at time 0). 

For any particular instance of the problem, finding a p(x,y,t) 
satisfying (3) and the boundary conditions is not difficult using numerical 
methods. Such procedures are routinely used in heat transfer problems 
to solve the diffusion equation (called the conduction equation or the 
Fourrier equation in the physics and engineering literature) to determine 
the temperature distribution across imperfectly conducting solids. In 
fact, Pitts and SissomI 1 977] give the example of a heated pipe in a 
square block of insulating material as one where the isotherms can be 
accurately estimated by hand plotting. 



RLthough the problem is not hard to solve numerically, the square 
boundary of R combined with a circular sensor tend to make the 
analytical solution difficult to obtain. Rnd without an analytical 
solution, it may be impossible to establish Sislioglu's observations in 
general. Making a change to the geometry, however, allows an 
analytical solution. Specifically, if the search region R is assumed to be 
a disk of area R instead of a square, then an exact solution in polar 
coordinates is possible. 

It is noted that the ray solution method described by Mangel[1981 ] 
could presumably be used to solve the diffusion equation, at least 
approximately, for a square search area. Such a solution was not 
attempted since an exact solution was available for the circular search 
area case. 



3. The Solution 

The disk-within-a-disk geometry has a radial symmetry, thus reducing 
the problem to one dimension, r. The new problem is to find a function 
p(r,t] satisfying 

(D/2) (d 2 p/dr 2 + (1/r)(3p/3r)J = 3p/3t, 

subject to 

9p/dr| r=Rft = 0, 
p(r,t) = 0, r4 R, 
p(r,0] = 1/fi, r>R, ‘ 

where R fl = (R/jr) 1/2 is the radius of the transformed (i.e., circular) search 
area fi. This problem has been solved in the physics literature. 



Muskat [ 1 93? ] arid Muskat[1934] la paper investigating the production 
rate of oil lueLLs) give the solution as 

00 

p(r,t) = _ (tT/A) Yj k n U(or n r) exp(-D or n 2 1/2) (9) 

n=1 

where 

k n = / (J 0 2 (or n R) ~ J j 2 ( Or n Rf|)), (10) 

U(or n r) = V^ornRft) Jo(» n r) - J 1 (af r ,R ft ) Y 0 (or n r), (11) 

and or n is the n th positive root of U(orR). (That is, the n* h smallest 
positive value of or satisfying U(orR) = 0.) Rlso, J 0 , J t , Y 0 , and Y, are 
respectively Bessel functions of the first kind order 0, first kind order 1, 
second kind order 0, and second kind order 1. 

The solution (9H11) does not appear particularly easy to evaluate or 
interpret, being an infinite sum of Bessel functions. Homever for large t 
the solution simplifies considerably. Rs t becomes large, (9) becomes 
Brk, /R) Ulc^r) exp(-D <*, 2 t/2), (12) 

where is the smallest positive value of or such that U(orR) = 0. The 
other exponential terms are ignored since they involve larger roots of 
U(orR). This means that as t becomes large, the decrease in p(r,t) for 
constant r becomes exponential at a rate of Da 2 / 2. 
lllhen p(r,t) is specified, PND(t) is then given by 

2* R fi 

J | p(r,t)rdrd8. (13) 



o R 



So for Large t, PND(t) becomes 



- (Zxr 2 k 1 /R) {Juicer) rdr} exp(-D <*, 2 t/2) 



(14) 



R 




= K exp (-D ar, 2 t/2), 

inhere K is a constant mhich depends on the problem geometry. 
Evaluation of the integral in (14) is straightforward given the change of 
variable u=o?, r and the identities 



Thus for large t, PND(t) decreases exponentially at the same rate as 
p(r,t). So Sislioglu's observation of exponentially decreasing PNDlt) (or, 
equivalently, a constant detection rate) appears asymptotically correct 
for large enough t. 

Using Muskat’s dimensionless notation, we can simplify the solution 
somewhat by defining 



Jx J 0 (x) dx = x J, (x) and Jx Y 0 (x) dx = x Y, (x). 



x, = or,R and p = Rf|/R. 

Then for large t, PND(t) can be written as 

K exp (-D X! 2 t/(2R 2 )), 
where x , is the first positive x satisfying 

Y,(xp) J 0 (x) - J,(xp) Y 0 (x) = 0; 

K is given by 



(16) 



-(277k, /(x, p 2 }) { J,(x, p) Y,(x,) - V,(x, p) J,(x,)}; 



(18) 



and k, is 



Uq(x,) J,(x, p) / (Jq 2 (x,) - J, 2 (x, p)). 



(19) 



In the next section, experimentally determined PND(t) from a Monte 
Carlo simulation of the diffusion process will be compared with (16H19), 
Sislioglu's approximation (1), and a slight modification of (1). 



4. Simulation Results and Conclusions 

Figure 1. is a plot of PND(t) determined by Monte Carlo simulation of 
the diffusion process for a square area R of size 10,000 square nautical 
miles (nm 2 ), a diffusion constant D of 100 nm 2 /hour, arid detection radii Fi 
varying from 28.21 nm (p=2) to 3.76 nm (p=15). The time increment for 
these simulations was 1 minute, which resulted in the x and y 
displacements of the target in each increment being selected from 
independent normal distributions of mean 0 and variance 100/60. 

Figure 1. illustrates the degree to which the decrease in PND(t) is 
exponential. Plotted on a log scale, PND(t) appears very nearly linear. 
For small t, however, the decrease in PND(t) is faster than exponential. 
This is seen more clearly in Figure 2, which is an enlargement of the 
upper left-hand corner of Figure 1. 

Sislioglu reported that PND(t) can be approximated by 

(1 - ttR 2 /R) exp(-24.7 RDt/fl 1 - 5 ). (20) 

The simulation results reported here indicate that a somewhat better 
approximation when p approaches 1 is 

(1 - 7TR 2 /R) exp(- 24.7 R0t/(R - TTR 2 ) 1 5 ). (21) 

That is, fl in (20) is replaced with fl - mrR 2 Figure 3. compares PN0(t) 
calculated by simulation with the estimates given by (21), (20), and the 
asymptotic estimate (16)-(19). The simulation data in Figure 3. are the 
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same as in Figures 1. and 2. except that only results for p of 2, 3, and 6 
are shown. In Figure 3. the simulated PND(t) is shown as a solid line, 
while the estimated PND(t) is shown as dashed. 

Figure 3.(0 indicates that, especially for small p, PND(t) determined 
by simulation does not decrease as rapidly as predicted by (16H19). The 
explanation is believed to be that for p sufficiently close to I, a circular- 
search region does not provide a reasonable approximation of a square 
region of the same area. To test this explanation, the three simulations 
plotted in Figure 3. were repeated with circular search areas. The 
results, shown in Figure 3.(D)., indicate a closer agreement between the 
theoretical and observed data. But still the fit is not exact. This is 
somewhat disturbing, but might be explained by the discrete manner in 
which the diffusion path is simulated. The simulated diffusion path is 
approximated by a series of points. Detection must occur exactly at the 
points and can not occur between them. The distance between adjacent 
points is the random variable (AX 2 +AY 2 ) 1/2 , where AX and AY are the 
independent, normally distributed x and y displacements. It is possible 
for the simulated path to jump across the edge of the detection disk 
without achieving detection, even though the line segment connecting 
the points lies partly on the disk. This will tend to reduce the simulated 
detection rate below that of the diffusion being approximated. 

Figure 4. shows a plot of x, vs. p for values of p from 2 to 15. By 
using (16), these values of Xj determine the theoretical asymptotic 
detection rate as 

D x 1 2 /(2R 2 ). (22) 

Table 1. lists the asympotic detection rates determined by (22), (21), (20), 




(y/ v y) oHa 



and an overall, (i.e., from time 0) best fit rate calculated by least- 
squares fitting of the simulation data. 



p 


X,(p) 


Dx 1 2 /(2R 2 ) 
(Eq. 22) 


24.7 RD/(R-JTR 2 ) 15 
(Eq. 21) 


24.7 RD/fl 1 - 5 
(Eq. 20) 


Simulation 
Best Fit 


2 


1.361 


.116 


.107 


.0697 


.102 


2.5 


.866 


.0736 


.0724 


.055? 


.0665 


3 


.626 


.0554 


.0554 


.0465 


.050? 


4 


.389 


.0381 


.0384 


.0384 


.0388 


6 


.218 


.0269 


.0242 


.0232 


.0250 


8 


.14? 


.021? 


.0178 


.0174 


.0202 


10 


.111 


.0194 


.0141 


.0139 


.01704 


15 


.0661 


.0154 


.00935 


.00929 


.0137 



Table 1. Detection Rates for Various p. 

The simulation data suggest the following two conclusions: 

a. Equations (16)-(19) give a reasonable estimate PND(t), but the fit 
deteriorates as p decreases to 1. 

b. For small p, (21) gives a better fit to the data than does (20). 
For large p, both (20) and (21) underestimate the observed asymptotic 
detection rate. 

LUe conclude with a few comments on random search. Rs mentioned 
earlier, Koopman's random search model predicts a detection rate of 
2Rv/R for searcher with speed v and detection range R conducting a 
random search of an area of size R. It seems reasonable that a 
diffusion path should be "random" enough for this model to be 
appropriate. In fact, we have seen that the detection rate, while not 



constant, approaches a constant value for large t. Setting 2Rv/R equal 
to the exponential term in (21) and solving for speed v gives 

12.35 RD/lR - 7TR 2 ) 15 

as an approximate equivalent speed for a searcher conducting a 
"diffusion search" of a square area. 
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